Convergence acceleration and stabilization for dynamical-mean-field-theory calculations 
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The convergence to the self-consistency in the dynamical-mean-field-theory (DMFT) calculations for models 
of correlated electron systems can be significantly accelerated by using an appropriate mixing of hybridization 
functions which are used as the input to the impurity solver. It is shown that the techniques and the past 
experience with the mixing of input charge densities in the density-functional-theory (DFT) calculations are 
also effective in DMFT. As an example, the increase of the computational requirements near the Mott metal- 
insulator transition in the Hubbard model due to critical slowing down can be strongly reduced by using the 
modified Broyden's method to numerically solve the non-linear self-consistency equation. Speed-up factors as 
high as 3 were observed in practical calculations even for this relatively well behaved problem. Furthermore, the 
convergence can be achieved in difficult cases where simple linear mixing is either not effective or even leads 
to divergence. Unstable and metastable solutions can also be obtained. We also determine the linear response 
of the system with respect to the variations of the hybridization function, which is related to the propagation of 
the information between the different energy scales during the iteration. 

PACS numbers: 71.27.+a, 71.30.+h. 72.15.Qm, 02.60.Cb 
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I. INTRODUCTION 

In many transition-metal, lanthanide, and actinide com- 
pounds the almost-localized d and / orbitals are par- 
tially filled and local magnetic moments are formed at low 
temperatures 1,2 . The competition between itinerancy and lo- 
cal electron-electron correlation effects gives rise to com- 
plex phase diagrams with different magnetic, charge-ordered, 
and superconducting phase s 1 ' 2 i 3 i 4 i 5 ' 6 . Simplified tight-binding 
models with short-range Coulomb interaction terms are com- 
monly used to study strong-correlation effects in such sys- 
tems. In the paradigmatic Hubbard mode l 7 ' 8 ' 9 i 1Q , the problem 
is reduced to a single-orbital description with purely on-site 
Coulomb repulsion. Hubbard-like models are thought to cor- 
rectly describe certain aspects of the itinerant ferromagnetism, 
the metal-insulator transitions 111 , and the high-temperature 
superconductivity^ 2 -. Despite intensive research, the prop- 
erties of the Hubbard model are not yet fully established. In 
the limit of infinite dimensions or high lattice connectivity the 
problem can be solved by the dynamical mean-field theory 
(DMFT) 'WVWWW 2 . In this limit, the self-energy 
becomes purely local and the bulk problem of correlated elec- 
trons maps exactly onto a quantum impurity model with a self- 
consistently-defined non-interacting bath of conduction elec- 
trons. In the DMFT, the spatial correlations are described in 
a mean-field way, however the local quantum fluctuations are 
taken into account exactly; as long as the effect under study 
is driven by local physics, the results of DMFT calculations 
are a good approximation for the properties of real (finite- 
dimensional) materials. 

In spite of the significant simplification of the full prob- 
lem within the DMFT, the solution of the effective quantum 
impurity problem is still challenging: it is by far the most 
computationally demanding part of the calculations. Among 
several impurity solvers in common use, the numerical renor- 
malization group (NRG) 23 ' 24 i 25 i 26 ! 27 i 28 i 29 i 30 is distinguished by 
its applicability to study the regime of very low tempera- 
tures directly in the thermodynamic limit. The convergence- 



acceleration approach proposed in the following is clearly ap- 
plicable to any impurity solver that may be used to solve the 
DMFT problem, however the discussion, the implementation 
details, and the test results are given for the NRG. Neverthe- 
less, the technique can easily be adapted for other solvers in a 
straight-forward manner. 

The input to a NRG calculation is the hybridization function 
r(w) which contains information about the density of states 
of the effective medium into which the impurity is embedded, 
while the output is an impurity spectral function A(u>), which 
is then used to compute the local lattice spectral function p(u>). 
The self-consistency is achieved when the two become equal, 
i.e. A(lu) = p(ui) within chosen accuracy, otherwise p(u>) 
is used to compute new hybridization function for the next 
DMFT iteration. In order to ensure the convergence, in some 
situations T(uj) from two consecutive iterations are linearly 
mixed to obtain the hybridization function which is used as the 
input to the NRG. A similar situation is well known in the field 
of quantum chemistry and electronic-structure calculations, in 
particular in the density-functional-theory (DFT) 31,32,33 where 
the quantity to be mixed is the charge density in space, n(r). 
In difficult cases (metallic surfaces, heterostructures, impu- 
rities in metals, systems near magnetic instabilities, etc.), 
the simple linear mixing procedure converges too slowly (or 
not at all), thus more sophisticated mixing approaches were 
devised 3 ^. In these schemes, a system of nonlinear equa- 
tions is solved iteratively by quasi-Newton-Raphson proce- 
dures or similar methoda 34 ' 35 ' 36 ! 3738 ' 39 . 40 . 41 ! 42 . 43 . 44 ! 45 . 46 . 47 . The 
solution of the system is equivalent to the Kohn-Sham vari- 
ational principle 48 . Such "advanced mixing" techniques are 
implemented with various degrees of sophistication in all DFT 
packages. They are stable, easy to implement and use, and 
they often have very high convergence rate. 

In this paper it is shown that the techniques and the past ex- 
perience with the advanced mixing schemes in DFT calcula- 
tions can also be applied to DMFT calculations. The advanced 
mixing greatly accelerates the convergence in many cases, for 
example near the Mott metal-insulator transition, where the 
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iteration converges very slowly due to critical slowing-down. 
It also ensures the convergence to unstable and metastable so- 
lutions, hence it can be applied to situations with multiple co- 
existing solutions. 

This work is structured as follows. In Sec. [IT] the DMFT 
self-consistency constraint is formulated as a sufficient condi- 
tion in the form of a system of equations. In Sec.|III]the modi- 
fied Broyden's iterative method for solving systems of nonlin- 
ear equations is briefly described, focusing on the implemen- 
tation with low storage requirements which is more suitable 
for large-scale problems 38 . In Sec. [TV] it is shown how the 
Broyden solver is incorporated into the DMFT loop and some 
further implementation details are given. In Sec.[Vjthe conver- 
gence properties of the linear and advanced mixing schemes 
are compared on the example of the Hubbard model for in- 
creasing electron-electron repulsion U. In Sec. [VTJ the Hub- 
bard model in external magnetic field is considered; in this 
case, the simple mixing is not always successful and the use 
of Broyden's method was found to be essential to obtain rapid 
convergence. In Sec. IVIII we study the response of the Hub- 
bard model with respect to small variations of the hybridiza- 
tion function; this response function is the equivalent of the 
Jacobian matrix of the system of self-consistency equations 
and describes the propagation of the information between var- 
ious energy scales during the DMFT iteration. Finally, in 
Sec. lVIIll some examples of the stabilization of otherwise un- 
stable fixed-point solutions are discussed. 



II. DMFT SELF-CONSISTENCY CONSTRAINT AS A 
SYSTEM OF NON-LINEAR EQUATIONS 
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(4) 



where ImA ff (w) = T a (uj) and the real part of A a (uj) can be 
obtained via the Kramers-Rronig transformation. In practice, 
the self-energy can be calculated more reliably and accurately 
as the ratio of two correlation functions 26 : the generalized 
Green's function F a (z) — ((d a n s ;d] T )) z over the Green's 
function G a (z), i.e. = UF a (ui)/G a (uj). The local 

lattice Green's function is 
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(5) 
(6) 
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where po (e) is the density of states (DOS) in the noninteract- 
ing model. The local lattice spectral function is then 



p a (u>) = Im [Gloc, CT (w + iS)} 



(8) 



The self-consistency condition 16 relates the local lattice 
Green's function Gi oc cr and the hybridization function T a as 

T <r (u) = -Tm[w-go t l(u)], (9) 
GoM") = G VoL +£<»• do) 



The single-orbital Hubbard model 7 ' 8 for electrons on a d- 
dimensional lattice 



(ij),<7 i,cr 



(with n i>(T = cl a c il<7 and = fj, - (cr/2)g(i B B) maps in 
the d — > oo limi t 16 ' 18 ' 49 onto the single-impurity Anderson 
models 



#SIAM = £d.an + Urn i n l + ^ [Vk^cl^da + H -' 

k,cr 



(2) 



with n„ = dld„ and 



n = n-f + nj_. The hybridization func- 



tion T a (u>) — J2k |Vfc,<r| 2 6(w — £fc.<r) contains full informa- 
tion about the coupling between the impurity and the effective 
non-interacting medium. From the calculated impurity spec- 
tral function 
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hn[G„{uj + iS)\ 



(3) 



where G a (z) = ((d a ; d' a )) z is the impurity Green's func- 
tion, one can extract the interaction self-energy T, a (uj) defined 



One DMFT cycle (which involves the numerical solution of 
-J/siam, the calculation of Gioc.o-, and the determination of the 
new hybridization function via Eq. (O) can be considered 
as a functional of the input hybridization function, i.e. 
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If the self-consistency has been established, the hybridization 
function is invariant (fixed point): 



rticw fT^old/, ,\~| j^old 

(7 \ (J \ / J (7 

Defining a mapping F as the difference 

F(r (T ) = r- w {r CT }-r ff , 



(12) 



(13) 



the approach to the self-consistency clearly corresponds to 
solving the system of equations 



F{T a ) = 0, 



(14) 



while a single DMFT step corresponds to applying F once 
to the hybridization function. Any solution of the equation 
dT4b is a possible physical state of the system since it satis- 
fies the self-consistency condition, albeit it is not necessar- 
ily the ground state: solutions corresponding to unstable and 
metastable states can also be found (see Sec. IVIIlt . It should 
be noted that Eq. ( fT4b is highly non-linear. 
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The usual DMFT iteration with no mixing corresponds to 
solving Eq. (fT4l i by a direct iteration, which often works since 
the mapping F behaves as a contraction in the vicinity of the 
solution (and often even far away from it). When F is not a 
contraction, however, this procedure will tend to diverge and 
more care is required. Usually it is sufficient to take an aver- 
age of two hybridization functions (the current and the previ- 
ous one): 

pinput,(m) _ a pncw,(m) -)- M _ Q ,)p in P ut >( m - 1 ) j (15) 

where a G [0 : 1] is the mixing parameter. It should be re- 
marked that this is fully analogous to simple charge mixing in 
the density-functional theory, where the charge density from 
the previous iteration n old (r) is admixed to the current one 
n new (r) as the true input to the next DFT iteration. Unfor- 
tunately, there are situations where this simple linear mixing 
approach fails even for small values of a. Furthermore, for 
very small a the approach to the self-consistency becomes 
prohibitively slow. In such situations, more sophisticated mix- 
ing approaches are required. In DFT, Broyden's method is 
commonly used. 



III. BROYDEN'S METHOD 

Let V be an TV-dimensional vector and F a mapping; the 
goal is to solve the system of equations F(V) = 0. The quasi- 
Newton-Raphson methods are iterative techniques in which 
the new approximation is given by 



y(m+l) _ y( m ) 



p(m) 



(16) 



where j( m ) is the Jacobian of the system at point V^™) and 
p(ro) _ p (\K m ->), xhe true Jacobian is unknown; a simple 
approximation is used for the initial Jacobian, for example a 
constant diagonal matrix 



-1, 



(17) 



which corresponds to simple linear mixing with a mixing pa- 
rameter a G [0 : 1]. The approximation is then improved 
by performing rank-1 updates as the iteration proceeds. It is 
more efficient to update directly the inverse of the Jacobian 

1 „ „35,37 



B ( m ) = - 

B (m+l) 

where 



^y(m) _ 
ApM _ 



y(m+l) _ y(m) 



|p(m+l) 
p(m+l) 



P(m)| 
p(m) 



|p(m+l) _ p(m) I 



AF (m) , 
(18) 

(19) 
(20) 



Vanderbilt and Louie have proposed a modified version of 
Broyden's method in which the information from all previ- 
ous iterations is incorporated when the current Jacobian is up- 
dated; this approach has better convergence properties and the 



Jacobian converges to the true Jacobian, which is not the case 
in the original Broyden's method which only uses the infor- 
mation from the most recent iteration to perform the update 36 . 
Srivastava has simplified the computational scheme so that 
only the input vectors V^™) and output vectors F^" 1 ) need to 
be stored, rather than the complete Jacobian matrix 37 . John- 
son combined the advantages of both schemes without any 
increase in complexity 38 . The final expressions for this modi- 
fied Broyden's method are as follows: 



m— 1 m— 1 



y(m+l) = y(m) + aF (m) « C i m '€ )uW 

(21) 



71=1 k=l 



with 



(m) = (aF^V F (m) 



U (n) = «AF (n) + AV (n) , 
and (to — 1) X (m — 1) dimensional matrices 



«( m ) _ 

Pk,n — 



A 



(m) 
k,n 



; w„(AF^) t AF^. 



(22) 
(23) 

(24) 
(25) 



Here 1 is a (to— 1) X (to— 1) dimensional identity matrix. The 
first two terms in Eq. (f2Tb correspond to simple linear mixing 
with parameter a, as described above, while the final term is 
a correction which takes into account the updates to the initial 
Jacobian. 

The weights w n (n = 1,2,...) are usually chosen to 
be equal to 1, while wo = 0.01322. For a suitable 
choice of weights, the modified Broyden's method becomes 
equivalent 39,45 to Pulay mixing scheme 42 or Anderson mixing 
scheme^ 4 -. 

The algorithm can be simply modified to use only a finite 
number of previous iterations to update the vector. This may 
be advantageous when the initial approximation for the vector 
is not very good. Alternatively, the Broyden mixing can be 
fully restarted after a given number of iterations. 



IV. INCORPORATION OF THE BROYDEN SOLVER INTO 
THE DMFT LOOP 

In the proposed convergence acceleration scheme for 
DMFT, the modified Broyden's method is used to refine the 
hybridization function T(ui) which is used as the input to the 
impurity solver. It should be remarked that this is not the only 
possibility: alternatively, one could also mix the self-energy 
£(u;). The choice depends somewhat on the problem and 
for numerical reasons one should in extreme cases preferably 
choose the quantity which is smoother as a function of the 
energy (for example, near the Mott transition on the metallic 
side the self-energy features sharp peaks while the hybridiza- 
tion function is rather smooth, whereas in the antiferromag- 
netic phase with small U the hybridization function contains 
sharp inverse-square-root singularities while the self-energy is 
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relatively smooth). In general, however, the two approaches 
are expected to be nearly equivalent. 

The Broyden solver is called just before the NRG, see 
Fig. 03 In a sense, the Broyden solver is effectively driving 
the DMFT loop in order to solve the equation 

F | r input : (m)| p(m+l) ^ -pinput,(m) 1 pinput.(m) q 

(26) 

see also Eqs. dT3l and 0141 . One cycle of the loop thus corre- 
sponds to applying once the mapping F to the hybridization 
function. 
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Figure 1 : (Color online) The DMFT loop using the numerical renor- 
malization group (NRG) as the impurity solver. The Broyden solver 
is incorporated in the loop as a correction step which modifies the 
input hybridization function in order to accelerate the convergence 
to the self-consistency. The new elements in the loop are shown in 
gray (red online). 

The vector V^" 1 ' corresponds to a discretized representa- 
tion of the continuous function T(ui). In the calculations pre- 
sented in this work, we used a geometric sequence of points 
u>±,n — i^mai™ with f2 max that exceeds the bandwidth of 
the noninteracting band by a factor of order 10 and I = 1.01; 
the dimension of the vectors was N = 3982 (and twice as 
large for spin-dependent problems). The same grid is used 
to sample all other functions, in particular the impurity spec- 
tral function A(ui) and the self-energy The Jacobian 
matrix in the Broyden solver was typically initialized with 
a = 1. The weights were chosen as wq = 0.01 and w n = 1 
for n > 1; setting wo to zero was found to have little ef- 
fect. It has been suggested that the weights w n be chosen as 
w n = (F( m '|F' m )) -1 / 2 , i.e. as the inverse root-mean-square 
difference of the function 38 . Numerical tests have shown that 
the improvement is only minor, if at all existing. 

During the initial steps it sometimes occurs that the result- 
ing r mput (oj) is not positive for all oj as the solver is overcom- 
pensating for the deviations. In such cases, the function was 
simply clipped to positive values. When the error vectors AF 
become smaller as the iteration proceeds, this is no longer a 
problem. The clipping performed during the initial iterations 



does not affect the final result. An alternative solution would 
be to revert to simple linear mixing in such instances. Yet an- 
other possible approach to enforce positivity of V would con- 
sist of working with In T instead. Unfortunately, this method 
was found to slow down the convergence significantly. 

It is necessary to store both r( m ) (oj) and r in P ut '( m ) (oj) for 
all A^stcps DMFT steps, thus the additional storage require- 
ments are of the order of N x iV s teps> which is not likely to 
pose difficulty. 

In calculations with fixed occupancy (rather than fixed 
chemical potential) it is important to store as an additional 
component of the vector V also the chemical potential [i that 
is being tuned. (This is actually true in general: all parameters 
varied in the iteration should appear in the Broyden process, 
so that the output of a single iteration is a smooth and uniquely 
defined function of the input vector V alone 40 .) In fact, the 
tuning of the parameter /i can be integrated in the Broyden 
solver with much fruition. 

In the NRG calculations performed for testing the method 
and presented in the following, the z-averaging 51 ' 52 ' 53 over 
N z = 8 values of the twist parameter was used in com- 
bination with an improved discretization scheme based on 
solving a differential equation to obtain the discretization 
coefficients^ 4 ^. The discretization parameter was A = 2, 
the truncation cutoff was set to E cuto s = IOljn (but no less 
than 500 and no more than 10000 states were used) and care 
was taken to truncate in a "gap" between clustered excitation 
levels. Spectral functions were computed using the density- 
matrix approach 56 and the self-energy trick 26 . Spectral infor- 
mation was extracted from both even and odd NRG iterations 
with a window parameter p — 2.3 54 . The broadening pro- 
cedure from Ref. 57 with a = 0.1 was used. The choice of 
NRG parameters appears to be important for the convergence: 
high-quality (smooth) results tend to be beneficial for the rate 
of convergence, while "rough" calculations sometimes lead 
to a stagnation of the convergence and oscillatory behavior. 
This is related to the assumption of differentiability of the 
mapping F. For the same reason, the calculations performed 
with larger broadening parameter will converge faster than 
high-energy -resolution calculations with much smaller broad- 
ening parameter. This is especially true when the hybridiza- 
tion function contains sharp features. 

The DMFT loop is terminated when two consecutive im- 
purity spectral functions A(oj) differ by no more than some 
chosen value: 

J |A (m) (u) - A (m - 1] {oj)\duj < A. (27) 

In practice it is found that this convergence test is more strin- 
gent when compared to an equivalent test for consecutive lo- 
cal lattice spectral functions p(uj), while comparing A(ut) and 
p(oj) at the same iteration gives absolute integrated errors 
somewhere between these two convergence tests. A typical 
convergence limit is A = 10 -6 . 

The stability of the converged solution can be tested by per- 
forming a few further DMFT iterations with the Broyden mix- 
ing turned off. From the solutions one can extract the domi- 
nant eigenvalue and eigenvector of the mapping F. This in- 
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formation is instrumental in assessing the physical stability of 
the solution and to determine the type of eventual instability. 
We return to these considerations in Sec. lVIIll 



V. ACCELERATION OF THE CONVERGENCE 

The acceleration of the convergence of the DMFT loop to- 
wards self-consistency was explored on the well-studied case 
of the Hubbard model at half-filling (/i = 0) in the param- 
agnetic regimo 16 ! 17 ! 18 ! 19 ! 22 ! 26 - 27 ! 29 ! 49 ! 58 . We study the Hubbard 
model on the Bethe lattice with infinite coordination number 
where 



(28) 



Here W is the width of the non-interacting conduction band. 
As the electron-electron repulsion U is increased, the char- 
acteristic three-peak structure emerges: two Hubbard bands 
and a quasiparticle peak at the Fermi level. As U approaches 
a critical value of U c /W ~ 1.46, the quasiparticle peak be- 
comes increasingly narrow and disappears 1 6 ' 26 ' 27,58 . Recent 
high-energy-resolution calculations have confirmed that the 
Hubbard bands have inner structure, in particular a peak at 
the inner edges 54,59 ' 60 ; this structure can be observed, for ex- 
ample, in the inset in Fig. [2] 




10 15 
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Figure 2: (Color online) Comparison of the convergence of the im- 
purity spectral function A{uS) in a calculation for the Hubbard model 
denned on the Bethe lattice in the paramagnetic phase at half-filling. 
We compare the simple mixing (here a = 1, i.e. the output from the 
previous iteration is used directly as the input for the new iteration, 
so this is actually direct iteration rather than mixing) and Broyden's 
mixing. The inset shows the converged density of states. 

In Fig. |2] we compare the approach to the self-consistency 
for the Hubbard model at fixed U/W = 1.2. The initial 
approximation for the local spectral function was the non- 
interacting DOS Po(ll>). Initially, both approaches are equiv- 
alent, since the starting approximation for the Jacobian is a 
diagonal matrix which corresponds to simple mixing. Since 
Po(w) is a rather crude approximation to the real density of 
states, Broyden's method is not expected to work much bet- 
ter than simple mixing for the first few steps; indeed, the 



errors are found to be even slightly higher. As can be seen 
in Fig. [3] the non-linear Broyden corrections are initially es- 
pecially large in the region of the emerging Hubbard bands, 
while at later iterations the most important contributions are to 
the inner-edge peaks in the Hubbard bands. Starting with iter- 
ation 9, the approximation to the self-consistent hybridization 
function becomes quite adequate, the updates to the Jacobian 
correspond to accurate refinements and the convergence accel- 
erates significantly. Both methods converge linearly, however 
the rate of convergence is much faster with Broyden's method. 
Superlinear convergence was never observed in practice. The 
linear rate of convergence in the example shown in Fig.|2]was 
/j w 0.15. When required, extremely good accuracy of the 
solution can thus be obtained with essentially no additional 
computational effort as compared to the direct iteration. 




Figure 3: The Broyden corrections for consecutive iterations; the 
curves are offset for clarity. The parameters are as in Fig. [2] 

We also determined the speed-up due to using the modi- 
fied Broyden's method as a function of the interaction strength 
U, Fig. |U As the Mott metal-insulator transition is ap- 
proached from below, the convergence becomes more difficult 
to achieve, which can be assigned to critical slowing down 
in the vicinity of quantum phase transitions 27 ' 61 ' 62 . Both ap- 
proaches are affected by this difficulty, however it is found 
that the relative speed-up in Broyden's method is an increas- 
ing function of U; for the range of parameter U considered 
in this work, the speed-up was up to a factor of 3 and it pre- 
sumably increases even further for U — ► U c . It should be re- 
marked that in these calculations it was possible to use a = 1 
(in other words, the direct DMFT iteration converges without 
any mixing), which is the most favorable situation. In prob- 
lems where mixing with small a is necessary, the speed-up 
factor is expected to be much higher. 

By performing calculation in the close vicinity of the Mott 
transition on the metallic side, we obtained an improved esti- 
mate of the critical value of U : 



U c /W = 1.459 . 



(29) 



It agrees very well with previous NRG calculations, where 
Uc/W = 1.47 was established 27 , and even better with 
the value obtained using projective self-consistent approach, 
U c /W = 1.46^. 
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Figure 4: (Color online) Comparison of the convergence as a func- 
tion of the interaction strength U/W. The convergence is defined to 
occur when two consecutive solutions for the density of states differ 
by no more than A = 1CP 6 (integrated absolute value of the differ- 
ence). The vertical dashed line corresponds to the point of the Mott 
metal-insulator transition at U c /W = 1.46 . 



VI. HUBBARD MODEL IN THE MAGNETIC FIELD 

The Hubbard model in a strong magnetic fiel d 63 ' 64 
has a metamagnetic response in a certain parameter 
regime: the magnetic susceptibility increases with the field 
strengt h 16 ' 63 ' 65 ' 66 ' 67 . The metamagnetic response is due to 
electron-electron interactions and, for sufficiently large U, it 
is driven mostly by field-induced quasiparticle mass enhance- 
ment (i.e. field-induced localization), however quasiparticle 
interactions also play a role 6 - 7 -. 

We consider the Hubbard model at half-filling and at zero 
temperature in a magnetic field. This problem is interesting 
for several reasons: 1) it is found that a DMFT iteration with 
simple mixing (and taking the noninteracting DOS as an ini- 
tial approximation) does not always converge in the presence 
of the magnetic field; 2) the structure at the inner-edge of the 
Hubbard band might be of magnetic origin, thus it can have 
non-trivial behavior in a finite magnetic field 54 - 5 ^; 3) the be- 
havior near the threshold to full polarization is not fully under- 
stood due to numerical difficulties in the transition regime 67 . 

The calculated spectral functions are presented in Fig. [5] 
The results agree with those shown in Ref. |67l however the 
energy resolution in our approach is sufficiently higher so that 
the inner structure in the Hubbard bands may be resolved 54 . 
As already established, when the magnetic field is increased 
the quasiparticle peak shifts away from the Fermi level and 
it narrows down, and the spectral weight is gradually trans- 
ferred to the lower Hubbard band of the majority spin 67 . With 
improved resolution, we can now also observe that the inter- 
nal structure of the Hubbard bands changes significantly with 
increasing field. When the magnetic field is increased past a 
transition value B c , a field-induced metal-insulator transition 
is induced 67 . 

For a system in the metallic regime, the Broyden's method 
converged rapidly, even when the non-interacting DOS was 
taken as the initial approximation, while linear mixing usu- 
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Figure 5: (Color online) Majority-spin spectral functions for the half- 
filled Hubbard model in a magnetic field. The arrows show the di- 
rection of the increasing magnetic field. 



ally led to oscillatory behavior. As in the Mott metal-insulator 
transition, the number of necessary iterations increases as the 
transition point is approached. On the insulating side, the con- 
vergence to the fully polarized solution was rapid for large 
fields, however the calculations in the vicinity of the tran- 
sition point were more difficult and it was necessary to ini- 
tialize the problem with the fully-polarized insulating spectral 
functions to ensure the convergence. The difficulties appear to 
stem from the fact that the non-interacting DOS for the Bethe 
lattice has square-root singularities at the band edges, while 
the Broyden method is premised on the differentiability of the 
mapping F. 

The inner structure in the Hubbard bands remains present 
even in the presence of the magnetic field, see Fig. [5] With 
increasing field, the lower Hubbard band of the majority-spin 
electrons becomes increasingly featureless and the inner-edge 
peak tends to disappear as we approach the transition to the in- 
sulating phase. The upper Hubbard band, however, appears to 
become more structured and distinctively asymmetric. Even 
at low fields there is some hint of further weak peaks within 
this band, which become more pronounced in the vicinity 
of the transition. In this regime, the electrons are already 
strongly polarized, thus majority-spin electrons in the upper 
Hubbard band cannot easily propagate since they reside on 
doubly-occupied sites surrounded predominantly by a ferro- 
magnetic background, thus their motion is strongly hindered 
by the Pauli exclusion principle and they become increasingly 
localized. This is to be contrasted with the holes in the lower 
Hubbard band which can easily propagate and do not feel the 
strong electron-electron interactions. 



VII. RESPONSE WITH RESPECT TO THE VARIATION 
OF THE INPUT HYBRIDIZATION FUNCTION 

Finding a good initial approximation to the Jacobian is not 
trivial, therefore a simple diagonal constant matrix is typically 
used, as in Eq. ( TP7I ). In band-structure calculations, the Jaco- 
bian is related to the dielectric tensor 34,37,48 , which makes it 
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possible to devise an improved initial approximation for the 
Jacobian based on the Thomas-Fermi screening theory (this 
procedure is called "preconditioning") 34 ! 37 ' 43 ! 46 . The Jacobian 
for the hybridization function in the DMFT loop is not related 
to some well-understood physical quantity in a simple way 
(see, however, the discussion of the Landau-Ginzburg func- 
tional Flg of the hybridization function discussed in Ref. |68| 
which is related to the self-consistency equation F(T) = 0). 
We may, however, study the properties of the Jacobian in the 
vicinity of the self-consistent solution r sc (o;) by performing 
calculations with slightly perturbed input hybridization func- 
tions: 



r in P ut ( w ) = r sc (w) + a - 



-[\u(u>/E)/b] 2 



(30) 



The perturbation takes the form of a log-Gaussian function 
centered at the energy E and of width b, similar to the com- 
monly used broadening kernel for producing smooth spec- 
tral functions in NRG (although the normalization factor 
differs) 29 . The weight a should be chosen small enough so 
that the response function 



R E (u>) = - [r out P ut (tj) - r sc (a 
a L 



(31) 



no longer depends on the value of a, but it must be large 
enough to prevent numerical artifacts. The width b should 
likewise be as small as possible, although its value is ulti- 
mately limited by the NRG broadening which restrains the 
energy resolution in r output (u;). The calculations were per- 
formed for a = 0.001 and b = 0.05. 
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Figure 6: (Color online) The response function Re(lo) for the Hub- 
bard model for the excitation energy E/W = 0.2 (graphically rep- 
resented as the vertical dashed line). Model parameters are the same 
as in Fig. [2] 

An example of the response function for the Hubbard model 
with intermediately strong interaction U/W — 1.2 is shown 
in Fig. [6] It reveals that a variation at a given excitation en- 
ergy E can lead to a complex response at all energies. (How- 
ever Re{u) vanishes in the uj — > limit 18 .) In simple linear 
mixing there is no exchange of information between different 
energies ("cross-talk"), thus it takes many DMFT iterations 



for reaching the self-consistency after a change has been im- 
posed. In the Broyden mixing, the application of the (approxi- 
mate) Jacobian effectively mixes the hybridization function at 
different energies, thereby accelerating the propagation of the 
information. 



VIII. UNSTABLE AND METASTABLE FIXED POINTS 

The concept of self-consistency is inseparably related to the 
concept of iteration; this is directly implied by the form of the 
self-consistency equation (fT2l . For this reason, the stability of 
the solutions (fixed points) is related to the eigenspectrum of 
the DMFT transformation, i.e. of the mapping F. Direct itera- 
tion can only be convergent if all the eigenvalues Aj of the lin- 
earization of F in the vicinity of the fixed point are strictly less 
than 1 in absolute value, while it will diverge when one or sev- 
eral eigenvalues are larger than one in absolute value, unless 
the solution space is constrained in such a way that the initial 
approximation for the solution has no components along the 
directions of the corresponding eigenvectors. For linear mix- 
ing with parameter a € [0 : 1] (note that a = 1 corresponds 
to direct iteration), the convergence criterium becomes 34 



|l-a(l-Ai)| < 1. 



(32) 



We denote by A max and A m i n the maximal and minimal eigen- 
value. If A max < 1 and A m j n > —1, the direct iteration with 
a = 1 will converge. If A max < 1 and A m ; n < — 1, a should 
be a < 2/(1 — A m j n ). Finally, if A max > 1 the inequality fl32"l > 
cannot be satisfied for any a E [0 : 1] and the linear mixing is 
of no help, thus the use of advanced mixing schemes becomes 
mandatory. 

As an example of a well-understood unstable solution, let 
us consider the instability of the paramagnetic solution of the 
Hubbard model at half-filling towards an antiferromagneti- 
cally ordered Neel ground stat o 18 ' 49 ' 69 ' 70 . Using Broyden's 
method, the paramagnetic (PM) solution can be stabilized in 
a calculation which in principle allows a symmetry broken 
state. The system drifts away, however, from the PM fixed 
point as soon as the Broyden mixing is turned off and eventu- 
ally it converges to an antiferromagnetic (AFM) solution, as 
illustrated in Fig. [7] The calculation was seeded with a pre- 
viously obtained self-consistent PM solution and iterated fur- 
ther without Broyden solver. The magnetization immediately 
starts to increase (Fig.|7p) and by the tenth iteration the spec- 
tral functions develop a narrow but sizeable singularity struc- 
ture (emerging spectral gap) in the quasiparticle peak, while 
the Hubbard bands start to become spin polarized (Fig. |7£). 
After 32 iterations, the result converged within 77 = 10 -6 to a 
stable self-consistent AFM solution shown in Fig.|7j5. 

It should be recalled that in a calculation where the SU(2) 
symmetry in spin space is explicitly maintained, the PM so- 
lution is stable and the mapping F is a contraction, as shown 
in Sec.|V] The eigenspectrum of the mapping F is not only a 
property of the physical model under consideration, but it also 
depends on the type of the long-range order allowed for in the 
DMFT equations, and to some degree even on the impurity 
solver used and on other details of the calculation (spectral 
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Figure 7: (Color online) Evolution of the paramagnetic solution for 
the half-filled Hubbard model after switching off the Broyden mix- 
ing, a) Initial paramagnetic (PM) spectral function, b) Resulting anti- 
ferromagnetic (AFM) spin-dependent spectral functions, c) Magne- 
tization and d) convergence as a function of the number of iterations, 
e) Difference between the spectral functions at iteration 9 and the ini- 
tial spectral function indicating the progressive breaking of the spin 
symmetry. 



broadening, discretization parameter, number of states kept, 
etc.). 

Another situation commonly encountered after switching 
off the Broyden mixing is the emergence of oscillatory solu- 
tions which never converge. This behavior can be observed 
for the Hubbard model at half-filling in a magnetic field (see 
Sec. I VP . The polarized fixed-point solution are found to be 
unstable and lead to oscillations between two almost fully 
spin-polarized (in the opposite directions) spectral functions, 
see Fig.[8h. The instability can be traced to a dominant eigen- 
value of A ~ —3.3 < — 1 (for the example in Fig. [8}, as ex- 
tracted from the ratio of differences between consecutive hy- 
bridization functions in the vicinity of the fixed point, Fig. ISJ). 
This solution could thus be stabilized using linear mixing with 
a < 0.23. It should be remarked that the solution was found to 
be unstable for all magnetic fields that yield a spin-polarized 
metallic solution, not only for weak fields where the system is 
known to be unstable toward the AFM solution. If the insta- 
bility of the fixed point is a true physical instability also for 
large magnetic fields, its nature is not very clear; it might cor- 
respond to canted ferromagnetism, a tendency towards forma- 
tion of spin density waves, or some other kind of incommen- 
surate order 71 . Since such states cannot be described by the 
formalism used, the iteration cannot converge. 

Since the fixed points T* of the self-consistency equa- 
tion F(T*) = are generally stationary points, rather than 
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Figure 8: Evolution of the solution for the Hubbard model at half- 
filling in an external magnetic field after switching off the Broyden 
mixing, a) Magnetization as a function of the number of iterations, b) 
ratio between the differences of consecutive hybridization functions, 
which provides an estimate for the dominant eigenvalue A ~ —3.3 
for B/W = 0.02. 



extrem a 68 ' 72 , it will be interesting to further clarify the re- 
lations between the stability of the DMFT iteration and the 
physical stability of the solution, as well as their relation to 
the eigenspectra of the mapping F in the vicinity of the so- 
lutions. As demonstrated, the proposed mixing technique can 
be a valuable tool for numerical studies of these questions, 
since it allows in principle to obtain all self-consistent solu- 
tions and (by turning the mixing off) to analyze the nature of 
their possible instabilities. 



IX. CONCLUSION 

It has been shown that the approach to the self-consistency 
can be greatly accelerated by reformulating the DMFT loop as 
an iterative method for solving a non-linear self-consistency 
equation using quasi-Newton-Raphson methods. The tests 
performed for the paradigmatic case of the Hubbard model 
at half-filling have shown that Broyden's method significantly 
outperforms simple linear mixing. The approach is fully gen- 
eral and it can be also applied when any other impurity solver 
(such as, for example, exact diagonalisation, DMRG, or quan- 
tum Monte Carlo) is used; it appears likely that similar speed- 
up factors could be achieved on equivalent problems. For 
particularly pathological situations, the improvement might 
be sufficient to bring previously forbidding problems within 
reach. This is particularly important near quantum phase tran- 
sitions, where reaching the convergence becomes problematic 
due to critical slowing down and the detailed behavior at the 
transition points is still a matter of controversy for many im- 
portant problems. The acceleration due to the use of Broy- 
den's method might be instrumental in answering some of 
these long-standing questions. In addition, the solver can be 
used to stabilize unstable fixed-point solutions and to study 
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their properties. Since the solver is robust, easy to implement 
and to incorporate in the DMFT cycle, there is little reason not 
to use it. 
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